
###############################################
########## JEPS OPTIMAL PERSUASION ############
###############################################

## Optimal Persuasion under Confirmation Bias: Theory and Evidence from a Registered Report
## Journal of Experimental Political Science
## Author: Love Christensen
## website: lovechristensen.com
## Contact: love.christensen@gmail.com

################################################
############ REPLICATION MATERIAL ##############
################################################

############### FILE 3: FIGURES ################

## It is important to run 1_experiment_prep_and_analysis.R before this file, since that file contains all the models which are called in creating the graphs

## GRAPH ORDER:

# 1. Figure 4: "Loess Curves of Updating"
# 2. Figure 5: "Non-Continuous Effects on Updating, Perceived Credibility and Support"
# 3. Appendix Figure 2: "Distribution of Prior Beliefs"
# 4. Appendix Figure 3: "Effect of Treatment on Prior Belief and Prior Certainty"
# 5. Appendix Figure 4: "Estimating the Share of Discounters"
# 6. Appendix Figure 5: "Effects of the Factorized Prediction Distance Variable"
# 7. Appendix Figure 6: "Copartisan Heterogeneity"
# 8. Appendix Figure 7: "Partisan Heterogeneity"
# 9. Appendix Figure 8: "Prior Certainty Heterogeneity"
# 10. Appendix Figure 10: "Treatment Assignment Distribution"


#################################################
##################### FIGURE 4 ##################
#################################################

## summarize data in bins
dat_bin <- dat %>%
  select(diff, preddistvalue) %>% 
  group_by(preddistvalue) %>%
  summarise(mean_diff = mean(diff, na.rm = TRUE), n = n()) %>% 
  arrange(preddistvalue)

## plot including prediction distance = 0
loess1 <- ggplot(dat, aes(preddistvalue, diff))  +
  geom_point(data = dat_bin, 
             mapping = aes(x = preddistvalue, y = mean_diff),
             alpha = 0.5) +
  geom_smooth(colour = "black")  +
  coord_cartesian(ylim = c(-1 , 0.5)) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        aspect.ratio=1,
        text = element_text(size=13))  +
  scale_y_continuous(breaks = seq(-1, 0.5, 0.5))  +
  xlab("Prediction Distance: $|x - \\mu_0|$ Millions of Jobs") +
  ylab("Updating ($\\mu_1 - \\mu_0$): Millions of Jobs") +
  ggtitle("a. Full Sample") +
  geom_hline(yintercept = 0, linetype = "dashed")

## plot exlcuding prediction distance = 0
loess2 <- ggplot(subset(dat, preddist != 0), aes(preddistvalue, diff))  +
  geom_point(data = dat_bin, 
             mapping = aes(x = preddistvalue, y = mean_diff),
             alpha = 0.5) +
  geom_smooth(colour = "black")  +
  coord_cartesian(ylim = c(-1 , 0.5)) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        aspect.ratio=1,
        text = element_text(size=13))  +
  scale_y_continuous(breaks = seq(-1, 0.5, 0.5))  +
  xlab("Prediction Distance: $|x - \\mu_0|$ Millions of Jobs") +
  ylab("Updating ($\\mu_1 - \\mu_0$): Millions of Jobs") +
  ggtitle("b. Excluding Null Prediction Distance") +
  geom_hline(yintercept = 0, linetype = "dashed")


fig.path <- paste('../output/figures/', 'fig_updating', '.tex', sep = "")
tikz(file = fig.path, width = 10, height = 5)
ggarrange(loess1, loess2, nrow = 1)
dev.off()

#################################################
##################### FIGURE 5 ##################
#################################################

## Extract coefficients and calculate confidence intervals
## I also reorder the cut variable so that everything is shown in the right order
## I use the same procedure for all coefficient plots

## updating
summary(m4_u)

dat_graph_u <- data.frame(summary(m4_u)$coefficients[-1, c("Estimate", "Std. Error")])
dat_graph_u <- dat_graph_u %>% 
  mutate(low95 = Estimate - Std..Error * 1.96,
         high95 = Estimate + Std..Error * 1.96,
         low90 = Estimate - Std..Error * 1.645,
         high90 =  Estimate + Std..Error * 1.645,
         value = c( "(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         value = factor(value, levels = value))

graph_u <- ggplot(dat_graph_u) +
          geom_linerange(aes(x = factor(value), 
                             ymin = low90,
                             ymax = high90),
                         size = 1.5) +
          geom_pointrange(aes(x = factor(value), 
                              y = Estimate, 
                              ymin = low95,
                              ymax = high95),
                          shape = 21,
                          fill = "white",
                          fatten = 6) +
          geom_hline(yintercept = 0, linetype = "dashed") +
          theme_minimal() +
          ggtitle("") +
          xlab("Prediction Distance Intervals") +
          ylab("Updating") +
  theme(aspect.ratio=1) +
  scale_y_continuous(breaks = seq(0.2, -0.8, -0.2),
                     limits = c(-0.75, 0.25))


## credibility
summary(m3_c)

dat_graph_c <- data.frame(summary(m3_c)$coefficients[-1, c("Estimate", "Std. Error")])
dat_graph_c <- dat_graph_c %>% 
  mutate(low95 = Estimate - Std..Error * 1.96,
         high95 = Estimate + Std..Error * 1.96,
         low90 = Estimate - Std..Error * 1.645,
         high90 =  Estimate + Std..Error * 1.645,
         value = c( "(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         value = factor(value, levels = value))

graph_c <- ggplot(dat_graph_c) +
  geom_linerange(aes(x = factor(value), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5) +
  geom_pointrange(aes(x = factor(value), 
                      y = Estimate, 
                      ymin = low95,
                      ymax = high95),
                  shape = 21,
                  fill = "white",
                  fatten = 6) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Credibility") +
  theme(aspect.ratio=1) +
  scale_y_continuous(breaks = seq(0.2, -0.6, -0.2),
                     limits = c(-0.75, 0.25))


## attitudes
summary(m1_a)

dat_graph_a <- data.frame(summary(m1_a)$coefficients[-1, c("Estimate", "Std. Error")])
dat_graph_a <- dat_graph_a %>% 
  mutate(low95 = Estimate - Std..Error * 1.96,
         high95 = Estimate + Std..Error * 1.96,
         low90 = Estimate - Std..Error * 1.645,
         high90 =  Estimate + Std..Error * 1.645,
         value = c( "(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         value = factor(value, levels = value))

graph_a <- ggplot(dat_graph_a) +
  geom_linerange(aes(x = factor(value), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5) +
  geom_pointrange(aes(x = factor(value), 
                      y = Estimate, 
                      ymin = low95,
                      ymax = high95),
                  shape = 21,
                  fill = "white",
                  fatten = 6) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Support") +
  theme(aspect.ratio=1) +
  scale_y_continuous(breaks = seq(0.2, -0.6, -0.2),
                     limits = c(-0.75, 0.25))

## combine into plot and write as tikz file
fig.path <- paste('../output/figures/', 'fig_outcomes', '.tex', sep = "")
tikz(file = fig.path, width = 9, height = 4)
ggarrange(graph_u, graph_c, graph_a, nrow = 1)
dev.off()

###########################################
############## APPENDIX FIGURE 2 ##########
###########################################

## "Distribution of Prior Beliefs"

prior_exp <- ggplot(dat, aes(tpp_prior_quant_1)) +
  geom_histogram(color="black", fill="white") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        aspect.ratio=1,
        text = element_text(size=13)) +
  scale_x_continuous(limits = c(-12.3, 12.3)) +
  labs(title="", x="Manufacturing Jobs Change", y = "Count")

prior_cert <- ggplot(dat, aes(tpp_prior_cert)) +
  geom_histogram(color="black", fill="white") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        aspect.ratio=1,
        text = element_text(size=13)) +
  scale_x_continuous(breaks=c(1:7)) +
  labs(title="", x="Certainty", y = "Count")


fig.path <- paste('../output/figures/', 'fig_priors', '.tex', sep = "")
tikz(file = fig.path, width = 9, height = 4)
ggarrange(prior_exp, prior_cert, nrow = 1)
dev.off()


###############################################
############### APPENDIX FIGURE 3 #############
###############################################

## "Effect of Treatment on Prior Belief and Prior Certainty"

## updating
dat_graph_priors <- data.frame(summary(m1_priors)$coefficients[-1, c("Estimate", "Std. Error")])

dat_graph_priors <- dat_graph_priors %>% 
  mutate(low95 = Estimate - Std..Error * 1.96,
         high95 = Estimate + Std..Error * 1.96,
         low90 = Estimate - Std..Error * 1.645,
         high90 =  Estimate + Std..Error * 1.645,
         preddist = seq(from = 1, to = 30, by = 1))

graph_priors <- ggplot(dat_graph_priors) +
  geom_linerange(aes(x = preddist, 
                     ymin = low90,
                     ymax = high90),
                 size = 1.2) +
  geom_linerange(aes(x = preddist, 
                     ymin = low95,
                     ymax = high95)) +
  geom_point(aes(x = preddist, 
                 y = Estimate),
             shape = 21,
             size = 2,
             fill = "white") +
  geom_smooth(data = dat_graph_priors, mapping = aes(preddist, Estimate),
              color = "gray") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance") +
  ylab("Prior Belief Expectation") +
  scale_y_continuous(breaks = seq(-2, 2, 0.5),
                     limits = c(-2, 2)) +
  scale_x_continuous(breaks = seq(0, 30, 5))

## credibility

dat_graph_priors2 <- data.frame(summary(m2_priors)$coefficients[-1, c("Estimate", "Std. Error")])

dat_graph_priors2 <- dat_graph_priors2 %>% 
  mutate(low95 = Estimate - Std..Error * 1.96,
         high95 = Estimate + Std..Error * 1.96,
         low90 = Estimate - Std..Error * 1.645,
         high90 =  Estimate + Std..Error * 1.645,
         preddist = seq(from = 1, to = 30, by = 1))

graph_priors2 <- ggplot(dat_graph_priors2) +
  geom_linerange(aes(x = preddist, 
                     ymin = low90,
                     ymax = high90),
                 size = 1.2) +
  geom_linerange(aes(x = preddist, 
                     ymin = low95,
                     ymax = high95)) +
  geom_point(aes(x = preddist, 
                 y = Estimate),
             shape = 21,
             size = 2,
             fill = "white") +
  geom_smooth(data = dat_graph_priors, mapping = aes(preddist, Estimate),
              color = "gray") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance") +
  ylab("Prior Belief Certainty") +
  scale_y_continuous(breaks = seq(-2, 2, 0.5),
                     limits = c(-2, 2)) +
  scale_x_continuous(breaks = seq(0, 30, 5))

fig.path <- paste('../output/figures/', 'fig_priors_treat', '.tex', sep = "")
tikz(file = fig.path, width = 8, height = 4)
ggarrange(graph_priors, graph_priors2, nrow = 1)
dev.off()


############################################
########### APPENDIX FIGURE 4 ##############
############################################

## "Estimating the Share of Discounters"

# Define the function for calculating what the share of discounters given the number of respondents, the initial slope and the end slope
# The process is described in Section C.3 of the appendix
f1 <- function(N, initial_slope, end_slope){
  n_exp_min <- 1
  n_exp_max <- N
  results_list <- list()
  for(i in n_exp_min:n_exp_max){
    x <- (end_slope - initial_slope * (4000 - i) / (4000) ) * (4000/i)
    results_list[[i]] <- data.frame(marg_disc = x, n_exp = i)
  }
  return(rbindlist(results_list) %>% mutate(n_exp_share = n_exp/4000))
}

## determine discounting for 4000 respondents, initial slope -0.57 and end slope 0.36
bounds_data <- f1(4000, -0.57, 0.36)


plot_bounds <- ggplot(data = bounds_data,
                      aes(x = n_exp_share, y = marg_disc)) +
  geom_line() +
  geom_hline(yintercept = 0.36, linetype = "dashed") +
  geom_hline(yintercept = 1.16, linetype = "dotted") +
  theme_minimal() +
  xlab("Share Discounting Respondents") +
  ylab("Marginal Persuasion of Discounting Respondents") +
  ggtitle("") +
  theme(aspect.ratio=1,
        plot.title = element_text(hjust = 0.5), 
        text = element_text(size=10)) +
  coord_cartesian(ylim = c(0, 5)) +
  scale_y_continuous(breaks = c(0, 0.36, 1, 1.16, 2, 3, 4, 5)) +
  theme(panel.grid.minor = element_blank())

fig.path <- paste('../output/figures/', 'fig_bounds', '.tex', sep = "")
tikz(file = fig.path, width = 3.75, height = 3.75)
plot_bounds
dev.off()


##############################################
############## APPENDIX FIGURE 5 #############
##############################################

## updating
dat_graph_u2 <- data.frame(summary(m10_u)$coefficients[-1, c("Estimate", "Std. Error")])
dat_graph_u2 <- dat_graph_u2 %>% 
  mutate(low95 = Estimate - Std..Error * 1.96,
         high95 = Estimate + Std..Error * 1.96,
         low90 = Estimate - Std..Error * 1.645,
         high90 =  Estimate + Std..Error * 1.645,
         preddist = seq(from = 1, to = 30, by = 1))

graph_u2 <- ggplot(dat_graph_u2) +
  geom_linerange(aes(x = preddist, 
                     ymin = low90,
                     ymax = high90),
                 size = 1.2) +
  geom_linerange(aes(x = preddist,
                     ymin = low95,
                     ymax = high95)) +
  geom_point(aes(x = preddist, 
                 y = Estimate),
             shape = 21,
             size = 2,
             fill = "white") +
  geom_smooth(data = dat_graph_u2, mapping = aes(preddist, Estimate),
              color = "gray") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance") +
  ylab("Updating (Change in Millions of Jobs)") +
  scale_y_continuous(breaks = seq(-2, 2, 0.5),
                     limits = c(-2, 2)) +
  scale_x_continuous(breaks = seq(0, 30, 5))

## credibility
dat_graph_c2 <- data.frame(summary(m10_c)$coefficients[-1, c("Estimate", "Std. Error")])
dat_graph_c2 <- dat_graph_c2 %>% 
  mutate(low95 = Estimate - Std..Error * 1.96,
         high95 = Estimate + Std..Error * 1.96,
         low90 = Estimate - Std..Error * 1.645,
         high90 =  Estimate + Std..Error * 1.645,
         preddist = seq(from = 1, to = 30, by = 1))

graph_c2 <- ggplot(dat_graph_c2) +
  geom_linerange(aes(x = preddist, 
                     ymin = low90,
                     ymax = high90),
                 size = 1.2) +
  geom_linerange(aes(x = preddist, 
                     ymin = low95,
                     ymax = high95)) +
  geom_point(aes(x = preddist, 
                 y = Estimate),
             shape = 21,
             size = 2,
             fill = "white") +
  geom_smooth(data = dat_graph_c2, mapping = aes(preddist, Estimate),
              color = "gray") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("")  +
  xlab("Prediction Distance") +
  ylab("Perceived Credibility of Prediction") +
  theme(aspect.ratio=1)  +
  scale_y_continuous(breaks = seq(-2, 2, 0.5),
                     limits = c(-2, 2)) +
  scale_x_continuous(breaks = seq(0, 30, 5))

## attitudes
dat_graph_a2 <- data.frame(summary(m10_a)$coefficients[-1, c("Estimate", "Std. Error")])
dat_graph_a2 <- dat_graph_u2 %>% 
  mutate(low95 = Estimate - Std..Error * 1.96,
         high95 = Estimate + Std..Error * 1.96,
         low90 = Estimate - Std..Error * 1.645,
         high90 =  Estimate + Std..Error * 1.645,
         preddist = seq(from = 1, to = 30, by = 1))

graph_a2 <- ggplot(dat_graph_a2) +
  geom_linerange(aes(x = preddist, 
                     ymin = low90,
                     ymax = high90),
                 size = 1.2) +
  geom_linerange(aes(x = preddist,
                     ymin = low95,
                     ymax = high95)) +
  geom_point(aes(x = preddist, 
                 y = Estimate),
             shape = 21,
             size = 2,
             fill = "white") +
  geom_smooth(data = dat_graph_a2, mapping = aes(preddist, Estimate),
              color = "gray") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance") +
  ylab("Support") +
  theme(aspect.ratio=1) +
  scale_y_continuous(breaks = seq(-2, 2, 0.5),
                     limits = c(-2, 2)) +
  scale_x_continuous(breaks = seq(0, 30, 5))

## combine into one plot
fig.path <- paste('../output/figures/', 'fig_outcomes_manycuts', '.tex', sep = "")
tikz(file = fig.path, width = 10, height = 3.5)
ggarrange(graph_u2, graph_c2, graph_a2, nrow = 1)
dev.off()

###############################################
################## FIGURE 6 ###################
###############################################

## "Copartisan Heterogeneity"

# The coefficient plots showing interaction effects are slightly more involved.
# This is because we want to plot the combined coefficients and not only the interaction coefficient
# To do this, we need the SE of the first and second variable and the covariance between them
# The function cov_finder takes as inputs the name of the two variables, and returns the covariance between these variables from the variance-covariance matrix
# All plots with interaction effects thereafter follow the same procedure

cov_finder <- function(namevar1, namevar2, vcov){
  vcov_element <- vcov[namevar1, namevar2]
  if(is.matrix(vcov_element)){
    vcov_element <- vcov_element[1,]
  }
  return(vcov_element)
}


## updating
summary(m1_up)

## set up first df for non-interaction coefficients
dat_pu1 <- data.frame(summary(m1_up)$coefficients[2:6, 1:2])
dat_pu1 <- dat_pu1 %>%
  mutate(varname1 = rownames(dat_pu1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

## set up second df for interaction coefficients
dat_pu2 <- data.frame(summary(m1_up)$coefficients[11:15, 1:2])
dat_pu2 <- dat_pu2 %>%
  mutate(varname2 = rownames(dat_pu2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(m1_up))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")

## merge the first df to the second df
## need estimates from both data frames to compute sum of coefficients
dat_pu2 <- left_join(dat_pu2, dat_pu1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se)) # SE for sum of coefficients

## calculate confidence intervals
dat_pu2 <- dat_pu2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         copartisan = "Yes")

## prepare dfs for appending
dat_pu2_app <- dat_pu2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         copartisan) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_pu1 <- dat_pu1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         copartisan = "No")

dat_pu1_app <- dat_pu1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         copartisan) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))
  
# append df1 and df2
dat_pu <- as.data.frame(rbind(dat_pu1_app, dat_pu2_app))
  
# plot them
graph_pu <- ggplot(dat_pu, aes(colour = factor(copartisan))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = copartisan),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Updating (Change in Millions of Jobs)") +
  theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Copartisan Sender ")  +
  scale_y_continuous(breaks = seq(-1.5, 0.5, 0.5),
                     limits = c(-1.5, 0.5))
  


## credibility
summary(m1_pc)

## set up dfs
dat_pc1 <- data.frame(summary(m1_pc)$coefficients[2:6, 1:2])
dat_pc1 <- dat_pc1 %>%
  mutate(varname1 = rownames(dat_pc1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_pc2 <- data.frame(summary(m1_pc)$coefficients[11:15, 1:2])
dat_pc2 <- dat_pc2 %>%
  mutate(varname2 = rownames(dat_pc2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(m1_pc))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")
dat_pc2 <- left_join(dat_pc2, dat_pc1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_pc2 <- dat_pc2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         copartisan = "Yes")

dat_pc2_app <- dat_pc2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         copartisan) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_pc1 <- dat_pc1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         copartisan = "No")

dat_pc1_app <- dat_pc1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         copartisan) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_pc <- as.data.frame(rbind(dat_pc1_app, dat_pc2_app))


graph_pc <- ggplot(dat_pc, aes(colour = factor(copartisan))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = copartisan),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Perceived Credibility") +
  theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Copartisan Sender ") +
  scale_y_continuous(breaks = seq(-1.5, 0.5, 0.5),
                     limits = c(-1.5, 0.5))

## attitudes
summary(m1_pa)

## set up dfs
dat_pa1 <- data.frame(summary(m1_pa)$coefficients[2:6, 1:2])
dat_pa1 <- dat_pa1 %>%
  mutate(varname1 = rownames(dat_pa1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_pa2 <- data.frame(summary(m1_pa)$coefficients[11:15, 1:2])
dat_pa2 <- dat_pa2 %>%
  mutate(varname2 = rownames(dat_pa2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(m1_pa))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")

dat_pa2 <- left_join(dat_pa2, dat_pa1, "varname1") 
dat_pa2 <- dat_pa2 %>%  
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_pa2 <- dat_pa2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         copartisan = "Yes")

dat_pa2_app <- dat_pa2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         copartisan) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_pa1 <- dat_pa1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         copartisan = "No")

dat_pa1_app <- dat_pa1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         copartisan) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_pa <- as.data.frame(rbind(dat_pa1_app, dat_pa2_app))


graph_pa <- ggplot(dat_pa, aes(colour = factor(copartisan))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = copartisan),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Support") +
  theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Copartisan Sender ") +
  scale_y_continuous(breaks = seq(-1.5, 0.5, 0.5),
                   limits = c(-1.5, 0.5))


## combine into one plot
fig.path <- paste('../output/figures/', 'fig_outcomes_partisan', '.tex', sep = "")
tikz(file = fig.path, width = 10, height = 5)
ggarrange(graph_pu, graph_pc, graph_pa, nrow = 1, common.legend = TRUE, legend="bottom")
dev.off()

##############################################
################# FIGURE 7 ###################
##############################################

## Partisan heterogeneity

## The principle is exactly the same as for the copartisan graph
## I first do the graphs for republicans, democrats and, lastly, inependents

############### REPUBLICANS ###################

################### UPDATING ##################

## model to graph
graph_mod <- m1_pr_u

## set up dfs
dat_temp1 <- data.frame(summary(graph_mod)$coefficients[2:6, 1:2])
dat_temp1 <- dat_temp1 %>%
  mutate(varname1 = rownames(dat_temp1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_temp2 <- data.frame(summary(graph_mod)$coefficients[8:12, 1:2])
dat_temp2 <- dat_temp2 %>%
  mutate(varname2 = rownames(dat_temp2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(graph_mod))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")
dat_temp2 <- left_join(dat_temp2, dat_temp1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_temp2 <- dat_temp2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         republicansender = "Yes")

dat_temp2_app <- dat_temp2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_temp1 <- dat_temp1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         republicansender = "No")

dat_temp1_app <- dat_temp1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_rep_upd <- as.data.frame(rbind(dat_temp1_app, dat_temp2_app))


graph_rep_upd <- ggplot(dat_rep_upd, aes(colour = factor(republicansender))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = republicansender),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Updating") +
  # theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Republican Sender ")+
  ggtitle("Republicans") +
  theme(plot.title = element_text(hjust = 0.5))  +
  scale_y_continuous(breaks = seq(-2, 1, 1),
                     limits = c(-2.5, 1.75))

################### CREDIBILITY ##################

## model to graph
graph_mod <- m1_pr_c

## set up dfs
dat_temp1 <- data.frame(summary(graph_mod)$coefficients[2:6, 1:2])
dat_temp1 <- dat_temp1 %>%
  mutate(varname1 = rownames(dat_temp1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_temp2 <- data.frame(summary(graph_mod)$coefficients[8:12, 1:2])
dat_temp2 <- dat_temp2 %>%
  mutate(varname2 = rownames(dat_temp2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(graph_mod))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")
dat_temp2 <- left_join(dat_temp2, dat_temp1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_temp2 <- dat_temp2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         republicansender = "Yes")

dat_temp2_app <- dat_temp2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_temp1 <- dat_temp1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         republicansender = "No")

dat_temp1_app <- dat_temp1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_rep_cre <- as.data.frame(rbind(dat_temp1_app, dat_temp2_app))


graph_rep_cre <- ggplot(dat_rep_cre, aes(colour = factor(republicansender))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = republicansender),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Credibility") +
  # theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Republican Sender ") +
  scale_y_continuous(breaks = seq(-2, 1, 1),
                     limits = c(-2.5, 1))


################### ATTITUDE ##################

## model to graph
graph_mod <- m1_pr_a

## set up dfs
dat_temp1 <- data.frame(summary(graph_mod)$coefficients[2:6, 1:2])
dat_temp1 <- dat_temp1 %>%
  mutate(varname1 = rownames(dat_temp1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_temp2 <- data.frame(summary(graph_mod)$coefficients[8:12, 1:2])
dat_temp2 <- dat_temp2 %>%
  mutate(varname2 = rownames(dat_temp2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(graph_mod))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")
dat_temp2 <- left_join(dat_temp2, dat_temp1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_temp2 <- dat_temp2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         republicansender = "Yes")

dat_temp2_app <- dat_temp2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_temp1 <- dat_temp1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         republicansender = "No")

dat_temp1_app <- dat_temp1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_rep_att <- as.data.frame(rbind(dat_temp1_app, dat_temp2_app))


graph_rep_att <- ggplot(dat_rep_att, aes(colour = factor(republicansender))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = republicansender),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Support") +
  #  theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Republican Sender ") +
  scale_y_continuous(breaks = seq(-1.5, 1, 0.5),
                     limits = c(-1.5, 1))


###############################################
################# INDEPENDENTS ################
###############################################

################### UPDATING ##################

## model to graph
graph_mod <- m1_pi_u

## set up dfs
dat_temp1 <- data.frame(summary(graph_mod)$coefficients[2:6, 1:2])
dat_temp1 <- dat_temp1 %>%
  mutate(varname1 = rownames(dat_temp1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_temp2 <- data.frame(summary(graph_mod)$coefficients[8:12, 1:2])
dat_temp2 <- dat_temp2 %>%
  mutate(varname2 = rownames(dat_temp2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(graph_mod))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")
dat_temp2 <- left_join(dat_temp2, dat_temp1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_temp2 <- dat_temp2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         republicansender = "Yes")

dat_temp2_app <- dat_temp2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_temp1 <- dat_temp1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         republicansender = "No")

dat_temp1_app <- dat_temp1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_ind_upd <- as.data.frame(rbind(dat_temp1_app, dat_temp2_app))


graph_ind_upd <- ggplot(dat_ind_upd, aes(colour = factor(republicansender))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = republicansender),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Updating") +
  #  theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Republican Sender ")+
  ggtitle("Independents") +
  theme(plot.title = element_text(hjust = 0.5))  +
  scale_y_continuous(breaks = seq(-2, 1, 1),
                     limits = c(-2.5, 1.75))

################### CREDIBILITY ##################

## model to graph
graph_mod <- m1_pi_c

## set up dfs
dat_temp1 <- data.frame(summary(graph_mod)$coefficients[2:6, 1:2])
dat_temp1 <- dat_temp1 %>%
  mutate(varname1 = rownames(dat_temp1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_temp2 <- data.frame(summary(graph_mod)$coefficients[8:12, 1:2])
dat_temp2 <- dat_temp2 %>%
  mutate(varname2 = rownames(dat_temp2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(graph_mod))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")
dat_temp2 <- left_join(dat_temp2, dat_temp1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_temp2 <- dat_temp2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         republicansender = "Yes")

dat_temp2_app <- dat_temp2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_temp1 <- dat_temp1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         republicansender = "No")

dat_temp1_app <- dat_temp1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_ind_cre <- as.data.frame(rbind(dat_temp1_app, dat_temp2_app))


graph_ind_cre <- ggplot(dat_ind_cre, aes(colour = factor(republicansender))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = republicansender),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Credibility") +
  # theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Republican Sender ")+
  scale_y_continuous(breaks = seq(-2, 1, 1),
                     limits = c(-2.5, 1))


################### ATTITUDE ##################

## model to graph
graph_mod <- m1_pi_a

## set up dfs
dat_temp1 <- data.frame(summary(graph_mod)$coefficients[2:6, 1:2])
dat_temp1 <- dat_temp1 %>%
  mutate(varname1 = rownames(dat_temp1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_temp2 <- data.frame(summary(graph_mod)$coefficients[8:12, 1:2])
dat_temp2 <- dat_temp2 %>%
  mutate(varname2 = rownames(dat_temp2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(graph_mod))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")
dat_temp2 <- left_join(dat_temp2, dat_temp1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_temp2 <- dat_temp2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         republicansender = "Yes")

dat_temp2_app <- dat_temp2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_temp1 <- dat_temp1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         republicansender = "No")

dat_temp1_app <- dat_temp1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_ind_att <- as.data.frame(rbind(dat_temp1_app, dat_temp2_app))


graph_ind_att <- ggplot(dat_ind_att, aes(colour = factor(republicansender))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = republicansender),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Support") +
  #  theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Republican Sender ") +
  scale_y_continuous(breaks = seq(-1.5, 1, 0.5),
                     limits = c(-1.5, 1))


###############################################
################### DEMOCRATS #################
###############################################

################### UPDATING ##################

## model to graph
graph_mod <- m1_pd_u

## set up dfs
dat_temp1 <- data.frame(summary(graph_mod)$coefficients[2:6, 1:2])
dat_temp1 <- dat_temp1 %>%
  mutate(varname1 = rownames(dat_temp1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_temp2 <- data.frame(summary(graph_mod)$coefficients[8:12, 1:2])
dat_temp2 <- dat_temp2 %>%
  mutate(varname2 = rownames(dat_temp2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(graph_mod))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")
dat_temp2 <- left_join(dat_temp2, dat_temp1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_temp2 <- dat_temp2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         republicansender = "Yes")

dat_temp2_app <- dat_temp2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_temp1 <- dat_temp1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         republicansender = "No")

dat_temp1_app <- dat_temp1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_dem_upd <- as.data.frame(rbind(dat_temp1_app, dat_temp2_app))


graph_dem_upd <- ggplot(dat_dem_upd, aes(colour = factor(republicansender))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = republicansender),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Updating") +
  #  theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Republican Sender ") +
  ggtitle("Democrats") +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(breaks = seq(-2, 1, 1),
                     limits = c(-2.5, 1.75))

################### CREDIBILITY ##################

## model to graph
graph_mod <- m1_pd_c

## set up dfs
dat_temp1 <- data.frame(summary(graph_mod)$coefficients[2:6, 1:2])
dat_temp1 <- dat_temp1 %>%
  mutate(varname1 = rownames(dat_temp1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_temp2 <- data.frame(summary(graph_mod)$coefficients[8:12, 1:2])
dat_temp2 <- dat_temp2 %>%
  mutate(varname2 = rownames(dat_temp2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(graph_mod))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")
dat_temp2 <- left_join(dat_temp2, dat_temp1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_temp2 <- dat_temp2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         republicansender = "Yes")

dat_temp2_app <- dat_temp2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_temp1 <- dat_temp1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         republicansender = "No")

dat_temp1_app <- dat_temp1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_dem_cre <- as.data.frame(rbind(dat_temp1_app, dat_temp2_app))


graph_dem_cre <- ggplot(dat_dem_cre, aes(colour = factor(republicansender))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = republicansender),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Credibility") +
  # theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Republican Sender ") +
  scale_y_continuous(breaks = seq(-2, 1, 1),
                     limits = c(-2.5, 1))


################### ATTITUDE ##################

## model to graph
graph_mod <- m1_pd_a

## set up dfs
dat_temp1 <- data.frame(summary(graph_mod)$coefficients[2:6, 1:2])
dat_temp1 <- dat_temp1 %>%
  mutate(varname1 = rownames(dat_temp1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_temp2 <- data.frame(summary(graph_mod)$coefficients[8:12, 1:2])
dat_temp2 <- dat_temp2 %>%
  mutate(varname2 = rownames(dat_temp2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(graph_mod))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")
dat_temp2 <- left_join(dat_temp2, dat_temp1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_temp2 <- dat_temp2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         republicansender = "Yes")

dat_temp2_app <- dat_temp2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_temp1 <- dat_temp1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         republicansender = "No")

dat_temp1_app <- dat_temp1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         republicansender) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_dem_att <- as.data.frame(rbind(dat_temp1_app, dat_temp2_app))


graph_dem_att <- ggplot(dat_dem_att, aes(colour = factor(republicansender))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = republicansender),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Support") +
  #theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Republican Sender ") +
  scale_y_continuous(breaks = seq(-1.5, 1, 0.5),
                     limits = c(-1.5, 1))



## Finally, we combine the nine plots into one large graph
graph_partisan_split <- ggarrange(graph_dem_upd, graph_ind_upd, graph_rep_upd,
                                  graph_dem_cre, graph_ind_cre, graph_rep_cre,
                                  graph_dem_att, graph_ind_att, graph_rep_att,
                                  common.legend = TRUE, legend="bottom", nrow = 3, ncol = 3)


## And save it as a tikz file
fig.path <- paste('../output/figures/', 'fig_partisan_split_outcomes', '.tex', sep = "")
tikz(file = fig.path, width = 12, height = 6)
graph_partisan_split
dev.off()


###############################################
############### APPENDIX FIGURE 8 #############
###############################################

## "Prior Certainty Heterogeneity"

## updating
summary(m1_cu)

## set up dfs
dat_cu1 <- data.frame(summary(m1_cu)$coefficients[2:6, 1:2])
dat_cu1 <- dat_cu1 %>%
  mutate(varname1 = rownames(dat_cu1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_cu2 <- data.frame(summary(m1_cu)$coefficients[11:15, 1:2])
dat_cu2 <- dat_cu2 %>%
  mutate(varname2 = rownames(dat_cu2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(m1_cu))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")
dat_cu2 <- left_join(dat_cu2, dat_cu1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_cu2 <- dat_cu2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         certain = "Yes")

dat_cu2_app <- dat_cu2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         certain) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_cu1 <- dat_cu1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         certain = "No")

dat_cu1_app <- dat_cu1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         certain) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_cu <- as.data.frame(rbind(dat_cu1_app, dat_cu2_app))


graph_cu <- ggplot(dat_cu, aes(colour = factor(certain))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = certain),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Updating (Change in Millions of Jobs)") +
  theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Certain") +
  scale_y_continuous(breaks = seq(-1.5, 0.5, 0.5),
                     limits = c(-1.5, 0.65))


## credibility
summary(m1_cc)

## set up dfs
dat_cc1 <- data.frame(summary(m1_cc)$coefficients[2:6, 1:2])
dat_cc1 <- dat_cc1 %>%
  mutate(varname1 = rownames(dat_cc1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_cc2 <- data.frame(summary(m1_cc)$coefficients[11:15, 1:2])
dat_cc2 <- dat_cc2 %>%
  mutate(varname2 = rownames(dat_cc2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(m1_cc))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")
dat_cc2 <- left_join(dat_cc2, dat_cc1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_cc2 <- dat_cc2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         certain = "Yes")

dat_cc2_app <- dat_cc2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         certain) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_cc1 <- dat_cc1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         certain = "No")

dat_cc1_app <- dat_cc1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         certain) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_cc <- as.data.frame(rbind(dat_cc1_app, dat_cc2_app))


graph_cc <- ggplot(dat_cc, aes(colour = factor(certain))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = certain),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Perceived Credibility") +
  theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Certain")  +
  scale_y_continuous(breaks = seq(-1.5, 0.5, 0.5),
                     limits = c(-1.5, 0.65))

## attitudes
summary(m1_ca)

## set up dfs
dat_ca1 <- data.frame(summary(m1_ca)$coefficients[2:6, 1:2])
dat_ca1 <- dat_ca1 %>%
  mutate(varname1 = rownames(dat_ca1)) %>%
  rename(estimate1 = "Estimate",
         se1 = "Std..Error")

dat_ca2 <- data.frame(summary(m1_ca)$coefficients[11:15, 1:2])
dat_ca2 <- dat_ca2 %>%
  mutate(varname2 = rownames(dat_ca2),
         varname1 = str_split(varname2, ":", simplify = TRUE)[, 1]) %>%
  mutate(cov_se = cov_finder(varname1, varname2, vcov(m1_ca))) %>%
  rename(estimate2 = "Estimate",
         se2 = "Std..Error")
dat_ca2 <- left_join(dat_ca2, dat_ca1, "varname1") %>%
  mutate(estimate_sum = estimate1 + estimate2,
         estimate_sum_se = sqrt(se1^2 + se2^2 + 2 * cov_se))

dat_ca2 <- dat_ca2 %>%
  mutate(low95 = estimate_sum - estimate_sum_se * 1.96,
         high95 = estimate_sum + estimate_sum_se * 1.96,
         low90 = estimate_sum - estimate_sum_se * 1.645,
         high90 = estimate_sum + estimate_sum_se * 1.645,
         certain = "Yes")

dat_ca2_app <- dat_ca2 %>%
  select(estimate_sum,
         low95,
         low90,
         high95,
         high90,
         varname1,
         certain) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate_sum) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_ca1 <- dat_ca1 %>%
  mutate(low95 = estimate1 - se1 * 1.96,
         high95 = estimate1 + se1 * 1.96,
         low90 = estimate1 - se1 * 1.645,
         high90 = estimate1 + se1 * 1.645,
         certain = "No")

dat_ca1_app <- dat_ca1 %>%
  select(estimate1,
         low95,
         low90,
         high95,
         high90,
         varname1,
         certain) %>%
  mutate(cut =  str_split(varname1, "\\)", simplify = TRUE)[, 2]) %>%
  rename(mid = estimate1) %>%
  select(-varname1) %>%
  mutate(cut = c("(5, 10]", "(10, 15]", "(15,20]", "(20, 25]", "(25, 30]"),
         cut = factor(cut, levels = cut))

dat_ca <- as.data.frame(rbind(dat_ca1_app, dat_ca2_app))


graph_ca <- ggplot(dat_ca, aes(colour = factor(certain))) +
  geom_linerange(aes(x = factor(cut), 
                     ymin = low90,
                     ymax = high90),
                 size = 1.5,
                 position = position_dodge(0.6)) +
  geom_pointrange(aes(x = factor(cut), 
                      y = mid, 
                      ymin = low95,
                      ymax = high95,
                      colour = certain),
                  shape = 21,
                  fill = "white",
                  fatten = 6,
                  position = position_dodge(0.6)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal() +
  ggtitle("") +
  xlab("Prediction Distance Intervals") +
  ylab("Support") +
  theme(aspect.ratio=1) +
  scale_colour_manual(values=c("black", "grey50"),
                      name = "Certain")  +
  scale_y_continuous(breaks = seq(-1.5, 0.5, 0.5),
                     limits = c(-1.5, 0.65))


## combine into one plot
fig.path <- paste('../output/figures/', 'fig_outcomes_certainty', '.tex', sep = "")
tikz(file = fig.path, width = 10, height = 5)
ggarrange(graph_cu, graph_cc, graph_ca, nrow = 1, common.legend = TRUE, legend="bottom")
dev.off()



#############################################
############ APPENDIX FIGURE 10 #############
#############################################

## "Treatment Assignment Distribution"

graph_assignment_to_treatment <- ggplot(subset(dat, !is.na(preddist)), aes(factor(preddist))) +
  geom_bar(color="black", fill="white") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        aspect.ratio=1,
        text = element_text(size= 9)) +
  scale_x_discrete(breaks = seq(0, 30, 5)) +
  labs(title="", x="Prediction Distance", y = "Count") +
  geom_hline(yintercept = mean(table(dat$preddist)), linetype = "dashed")


fig.path <- paste('../output/figures/', 'fig_treatassign', '.tex', sep = "")
tikz(file = fig.path, width = 4, height = 3)
graph_assignment_to_treatment
dev.off()
